Andrea Luciani is a Technical Advisor for the Directorate General for Economics, Statistics and Research at the Bank of Italy, and co-author of the bimets package.

In our first post we presented an approach on how to estimate and simulate Structural Equation Models (SEM) in R. In particular, we focused on R tools that allow users to forecast advanced simultaneous equations models having linear restrictions on coefficients, error autocorrelation on residuals, and conditional equation evaluation.

The simulation of these econometric models often requires using iterative algorithms in ways that are well beyond what the ts(), lm() and predict() functions were designed to do; therefore, in our first post, we made a forecasting exercise of a Klein model by using the deterministic simulation procedure available into the bimets package.

However, deterministic algorithms have significant limitations when applied in forecasting exercises. Forecasts produced by structural econometric models are subject to several sources of error, such as random disturbance term of each stochastic equation, errors in estimated coefficients, errors in forecasts of exogenous variables.

The forecast error depending on the structural disturbances can be analyzed by using a stochastic simulation approach.

In stochastic simulations, the structural disturbances are given values with specified stochastic properties. The error terms of the estimated equations are appropriately perturbed. Technical identities and exogenous variables can also be perturbed by disturbances with specified stochastic properties. Programmatically, a straightforward approach to solve the problem is a Monte-Carlo method: the model is solved for each data set with different values for the disturbances. Finally, forecast statistics (e.g. mean, standard deviation, etc.) can be computed for each simulated endogenous variable.

Using the Klein model defined in our first post, we can attempt to complete a US Gross National Product forecasting exercise, given a user-specified uncertainty on Consumption and Government Expenditure time series.

The Klein model will be perturbed by applying a normal disturbance to the endogenous Consumption behavioral cn in year 1942, and a uniform disturbance to the exogenous Government Expenditure time series g along all the forecast time range. The normal disturbance applied to the cn stochastic equation has a zero mean and a standard deviation equal to its regression standard error, thus roughly replicating the regression error during the current perturbation (not accounting for inter-equations cross-covariance).

#we want to perform a stochastic forecast of the GNP up to 1944
#we will add normal disturbances to endogenous Consumption 'cn' 
#in 1942 by using its regression standard error
#we will add uniform disturbance to exogenous Government Expenditure 'g'
#in whole TSRANGE

myStochStructure <- list(
  cn=list(
    TSRANGE=c(1942,1,1942,1),
    TYPE='NORM',
    PARS=c(0,kleinModel$behaviorals$cn$statistics$StandardErrorRegression)
  ),
  g=list(
    TSRANGE=TRUE,
    TYPE='UNIF',
    PARS=c(-1,1)
  )
)
 
#model stochastic forecast 
kleinModel <- STOCHSIMULATE(kleinModel
                            ,simType='FORECAST'
                            ,TSRANGE=c(1941,1,1944,1)
                            ,StochStructure=myStochStructure
                            ,StochSeed=123
                            ,quietly=TRUE
)

#print mean and standard deviation for the forecasted GNP
with(kleinModel$stochastic_simulation,TABIT(y$mean, y$sd))
## 
##       Date, Prd., y$mean         , y$sd           
## 
##       1941, 1   ,  125.5045      ,  4.250935      
##       1942, 1   ,  173.2946      ,  9.2632        
##       1943, 1   ,  185.9602      ,  11.87774      
##       1944, 1   ,  141.0807      ,  11.6973

Another limitation of the deterministic simulation procedure relies on the fact that analysts often do not want to perform a mere forecasting exercise. Usually, a performance measure (i.e. objective-function) is assigned to an econometric model, depending on the value of forecasted endogenous variables; thus, analysts try to enhance this measure by fine-tuning exogenous variables (e.g. the instruments) of the model (e.g. policy analysis).

Generally speaking, this is an optimal control exercise: the optimization consists of maximizing an objective-function, depending on (forecasted) endogenous variables, given a set of user-constraints plus the constraints imposed by the econometric model equations.

In the following exercise, we will use the bimets OPTIMIZE() procedure which allows R users to perform optimal control exercises on simultaneous equations models. The exercise will be completed by using the following code on the previously defined Klein model, and we will assume that:

  1. The objective-function definition is: \(f(y, cn, g) = (y-110)+(cn-90)*|cn-90|-\sqrt{g-20}\) given \(y\) as the simulated Gross National Product, \(cn\) as the simulated Consumption and \(g\) as the exogenous Government Expenditure. The basic idea is to maximize Consumption, and secondarily the Gross National Product, while reducing the Government Expenditure;

  2. The instrument variables (i.e. INSTRUMENT in following code snippet) are the \(cn\) Consumption “booster” (i.e. the add-factor, not to be confused with the simulated Consumption in the objective-function) and the \(g\) Government Expenditure, defined over the following domains: \(cn \in (-5,5)\), \(g \in (15,25)\);

  3. The following restrictions are applied to the INSTRUMENT: \(g + cn^2/2 < 27 \wedge g + cn > 17\), given again \(cn\) as the Consumption “booster” (i.e. the add-factor) and \(g\) as the Government Expenditure;

The scatter plot at the end of this post is populated with 2916 objective-function stochastic realizations; the 210.58 local maximum is stored in the kleinModel$optimize$optFunMax variable, while the instruments that allow local maximum to be achieved are stored in the kleinModel$optimize$INSTRUMENT variable. Code follows:

#we want to maximize the non-linear objective function:
#f()=(y-110)+(cn-90)*ABS(cn-90)-(g-20)^0.5
#in 1942 by using INSTRUMENT cn in range (-5,5) 
#(cn is endogenous so we use the add-factor)
#and g in range (15,25)
#we will also impose the following non-linear restriction:
#g+(cn^2)/2<27 & g+cn>17

#define INSTRUMENT and boundaries
myOptimizeBounds <- list(
    cn = list( TSRANGE = TRUE
            ,BOUNDS = c(-5,5)),
     g = list( TSRANGE = TRUE
            ,BOUNDS = c(15,25))
)

#define restrictions
myOptimizeRestrictions <- list(
    myRes1=list(
         TSRANGE = TRUE
        ,INEQUALITY = 'g+(cn^2)/2<27 & g+cn>17')
)

#define objective function
myOptimizeFunctions <- list(
    myFun1 = list(
         TSRANGE = TRUE
        ,FUNCTION = '(y-110)+(cn-90)*ABS(cn-90)-(g-20)^0.5')
)

#Monte-Carlo optimization by using 10000 stochastic realizations
#and 1E-4 convergence criterion 
kleinModel <- OPTIMIZE(kleinModel
                      ,simType = 'FORECAST'
                      ,TSRANGE=c(1942,1,1942,1)
                      ,simConvergence= 1E-4
                      ,simIterLimit  = 1000
                      ,StochReplica  = 10000
                      ,StochSeed = 123
                      ,OptimizeBounds = myOptimizeBounds
                      ,OptimizeRestrictions = myOptimizeRestrictions
                      ,OptimizeFunctions = myOptimizeFunctions
                      )
## OPTIMIZE(): optimization boundaries for the add-factor of endogenous variable "cn" are (-5,5) from year-period 1942-1 to 1942-1.
## OPTIMIZE(): optimization boundaries for the exogenous variable "g" are (15,25) from year-period 1942-1 to 1942-1.
## OPTIMIZE(): optimization restriction "myRes1" is active from year-period 1942-1 to 1942-1.
## OPTIMIZE(): optimization objective function "myFun1" is active from year-period 1942-1 to 1942-1.
## 
## 
Optimize:     100.00 %
## OPTIMIZE(): 2916 out of 10000 objective function realizations (29%) are finite and verify the provided restrictions.
## ...OPTIMIZE OK
#print local maximum
kleinModel$optimize$optFunMax
## [1] 210.5755
#print INSTRUMENT that allow local maximum to be achieved
kleinModel$optimize$INSTRUMENT
## $cn
## Time Series:
## Start = 1942 
## End = 1942 
## Frequency = 1 
## [1] 2.032203
## 
## $g
## Time Series:
## Start = 1942 
## End = 1942 
## Frequency = 1 
## [1] 24.89773

Disclaimer: The views and opinions expressed in this page are those of the author and do not necessarily reflect the official policy or position of the Bank of Italy. Examples of analysis performed within these pages are only examples. They should not be utilized in real-world analytic products as they are based only on very limited and dated open source information. Assumptions made within the analysis are not reflective of the position of the Bank of Italy.